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FINITE ELEMENT COMPUTER PROGRAM TO ANALYZE 
CRACKED ORTHOTROPIC SHEETS 


By 

Chorng-Shin Chu 
J. M. Anderson* 
W. J. Batdorf 
J. A. Aberson** 


SUMMARY 


The objective of this study was to develop a finite-element computer program that 
performs a two-dimensional elastostatic analysis of plane anisotropic homogeneous sheets 
with a through-the-th.ickness crack. The program includes special crack-tip elements that 
account for the singular stress fields associated with crack opening (Mode I) and sliding 
(Mode II) displacements at the crack tips. These special crack-tip elements provide a new 
tool for computing stress- intensity factors, and they can be used for predicting the crack 
propagation for a damaged structure. 

Two types of crack elements have been developed. One element has 8 nodes and is 
restricted to symmetric applications; the other element has 10 nodes and is capable of repre- 
senting the crack-tip neighborhood when both the crack opening and sliding modes of defor- 
mation occur. The 8-node symmetric cracked element takes full advantage of the symmetry, 
so that only one-half of a configuration which is symmetric about a line containing the 
crack need be modeled. However, the 10-node unsymmetric cracked element has much 
wider usage in practical fracture-mechanics applications. 

These cracked elements can exhibit either anisotropic or isotropic behavior. (Ortho- 
tropic materials are considered a special case of anisotropic material.) A stress function in 
a half-power series form appropriate to the crack tip neighborhood is chosen to formulate 
the stiffness matrix for the anisotropic cracked elements, and a Williams' series stress func- 
tion is used for the isotropic cracked elements. 


Consultant, Associate Professor, School of Engineering Science and Mechanics, Georgia 
Institute of Technology. 

**Consultant, Assistant Professor, School of Engineering Science and Mechanics, Georgia 
Institute of Technology. 



The computer program contains an axial element, a linear spring element, a triangular 
element, and isotropic and anisotropic cracked elements. Thermal-strain analysis capability 
is included. Extensive tests during the development stage have been made in order to vali- 
date the program against known solutions. To illustrate the capabilities of these two types 
of cracked elements, a number of selected sample problems are presented in this report for 
demonstrations. The computer program has proved to be very efficient. 


INTRODUCTION 


The finite-element method is one of the most effective approaches available for plane 
elasticity problems having irregular boundaries or discontinuous boundary conditions. Early 
efforts to bring this method to bear on crack problems depended on the use of many conven- 
tional elements around the crack tip in an attempt to represent the extreme stress gradients 
there. Stress- intensity factors were estimated from results obtained with such models either 
by extrapolating crack-opening displacement (ref. 1) or by numerically computing the vari- 
ation in strain energy with crack length (ref. 2). Such conventional methods have been 
found to be limited and uneconomical. For example, Oglesby and Lomacky (ref. 3) have 
indicated that the maximum permissible element size necessary to ensure acceptable accu- 
racy (5% error or less) is on the order of 1/500 of the crack half length. 

To circumvent this economic problem, development and research efforts have turned 
toward formulating elements which contain the crack-tip stress singularity. These special 
singularity elements, usually referred to as cracked finite elements, represent a significant 
improvement in both accuracy and economy in comparison with conventional methods. Many 
cracked finite elements developed to date (ref. 4-8) incorporate only the singular term in 
the series expansion for the crack-tip stress field. Admittedly, this term dominates all 
others near the crack tip, but to guarantee that the nonsingular contributions are comparably 
negligible, the neighborhood represented by the cracked element must be quite small, and 
the problem of economy arises again. Moreover, the neighboring conventional elements in 
such instances are drawn very near the crack tip again, and concern over their capability 
to represent the stress field adequately has led some investigators to introduce special 
"border elements" having a higher degree of sophistication than that routinely required for 
plane elasticity problems. 

Wilson (ref. 9) has developed a cracked finite element that makes use of the first four 
terms in the expansion for the crack-tip stress field. He reports accurate results when this 
element is used in conjunction with a fairly modest number of conventional elements. 
Wilson's element, however, has the disadvantage of being semi-circular and hence is some- 
what awkward to use with conventional elements, which almost always have straight bound- 
aries. Moreover, Wilson's element (as well as some others previously referenced) has fewer 
degrees of freedom than are needed for independence of the nodal displacements. This 
requires that the stiffness matrix of the cracked element receive special attention in form- 
ing the stiffness matrix of the assembly. 
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Early in 1973, the Lockheed-Georgia Company completed the development of isotropic 
versions of two high-order cracked finite elements, i.e., elements that incorporate many 
terms in the expansion for the crack-tip stress field. This feature permits very accurate esti- 
mates of stress-intensity factors with relatively coarse finite-element grids. Both high-order 
cracked elements, one for symmetric (Mode I) applications and one for unsymmetric (Mode I 
and Mode II) applications, have a perfect balance between actual degrees of freedom and 
the number of nodal displacement components. Thus, the numerical analyst adds the cracked 
element to an assembly in exactly the same way that he adds a conventional element. The 
shape of each element was chosen to make it fit conveniently in models making use of the 
most widely used membrane element: the constant-strain triangle. The remarkable accuracy 
obtained with both elements in isotropic applications using very coarse finite-element repre- 
sentations (ref. 10) would seem to indicate that the periodic concern of some investigators 
over the displacement incompatibility that exists at the interface of the high-order cracked 
element and a conventional element is largely academic. 

The continually increasing use of fiber-reinforced composites in aerospace applications, 
coupled with a growing confidence in the ability of linear elastic fracture mechanics (LEFM) 
to predict the growth rate and stability of cracks in isotropic materials, has lately resulted in 
considerable interest in the prospects of successfully applying LEFM to anisotropic materials 
(ref. 1 1). This report contains the analytical background and an account of several numeri- 
cal modifications and additions required to effect the following principal extensions of 
Lockheed-Georgia's existing analysis capability for cracked isotropic structures: 

o Axial element 

o Anisotropic constant-strain triangle 
o Anisotropic cracked elements 

o Automatic node sequencing to minimize bandwidth 
o Input-data generator 
o Thermal stress analysis capability 


SYMBOLS 


w 

l P r, ! 

0] 

HI 


Vector of forces for uncoupled structure 
Vector of displacements for uncoupled structure 

Vector of thermal forces for complete restraint of uncoupled structure 

Uncoupled stiffness matrix 

Vector of displacements for coupled structure 
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T r9 
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Vector for forces for coupled structure 
Compatibility matrix for displacements 
Coupled stiffness matrix 

Vector of thermal forces for complete restraint of coupled structure 
Stress vector 
Strain vector 

Elastic matrix for Hooke's law in local axis 
Elastic matrix for Hookes law in material axis 
Transformation matrix for elastic coefficients 
Vector representing thermal expansion coefficients 
Rectangular coordinates 
Polar coordinates 

Normal components of stress parallel to x-, y-, and z-axes 
Shearing-stress components in rectangular coordinates 
Unit elongations in x-, y- , and z-directions 
Shearing-strain components in rectangular coordinates 
Radial and tangential normal stresses in polar coordinates 
Shearing stress in polar coordinates 

Radial and tangential unit elongation in polar coordinates 

Modulus of elasticity in tension and compression 

Modulus of elasticity in shear 

Poisson's ratio 

Strain energy 

Temperature 
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T 


u 


Stress function 


a Half crack length 

L Total height of cracked tension plate 

W Total width of cracked tension plate 

w = W/2 

a/w Crack length aspept ratio 
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BASIC EQUATIONS FOR FINITE ELEMENT PROGRAM 


Matrix Equations for Stiffness Method 

The finite element analysis program is based on the direct stiffness method. With this 
approach, a complex structure is idealized as an assembly of simple elements; for example: 
axial elements, triangular membrane elements (isotropic and anisotropic), and 8- and 10- 
node cracked elements. Each of these elements, separately, can be analyzed without diffi- 
culty. By expressing the various mathematical relationships in matrix form, it is possible to 
assemble automatically a large number of discrete elements to simulate almost any complex 
structure. The basic matrix operations that are executed by the program are illustrated here. 
These equations include thermal effects. 

In matrix form, the uncoupled force-displacement equation for a single element or the 
entire structure can be written as 

{p}=[k]{ P } + {P rt } (i) 

where P and p are vectors of forces and displacements, respectively, and k is the uncoupled 
stiffness matrix. The term P r1 . is a vector of thermal forces equivalent to the forces produced 
when each element is completely restrained. 

From compatibility considerations we can define tb such that 


H-MH-Mii 


where q is a vector of node point displacements for the coupled structure. In other words 
the p vector is referred to the local uncoupled system and q to the coupled structure de- 
fined in some global reference frame. The matrix ij) consists of direction cosines and is 
determined by the topology of the structural model. The vector q can then be expressed in 
terms of unknown displacements q Q and known displacements q^. From the principal of 
virtual work it can be shown that 


{o}=[/k *]{,}+ |7]{p rt } 

(3) 

H- WM + k,} 

(4) 


Equation (4) is then the force-displacement equation for the coupled or global system, 
where Q and q are forces and displacements, respectively. The vector Q r j. represents 
forces for complete restraint for the global system. Equation (4) can be expanded as 
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where Q a and q a are applied forces and unknown displacements; and and q^ are 
reaction forces and known displacements. To determine the elastic displacements for a hot 
structure it is then necessary to solve the following set of linear equations: 

{ Q a- K ab%- Q ar,} = [ K a=]k} < 6 > 

Once the displacements are known, the loads in each discrete element can then be deter- 
mined from the following equation 

M41WM+W (7 ) 

The reaction forces, alonq with the equilibrium check, can be determined from Equation (8) 

M- MM « 

References (12) and (13) discuss the various concepts of matrix analysis of structures. 
The derivations of the element stiffness matrices, except for the cracked elements, are out- 
lined in the Appendix. 


Transformation of Elastic and Thermal Coefficients 


For a two-dimensional elastic plate the relationship between stress and strain can be 
written as 


M= !>.]{<} 


(9) 


For an isotropic material and plane-stress conditions. 




0 1 



0 


L 0 0 G J 


( 10 ) 


In the general case of anisotropic materials the only requirement is that the matrix of co- 
efficients be symmetric, that is 


OD 


M - 


AAA 
A 11 12 16 

A 12 A 22 A 26 

L A 16 A 26 A 66 


The six coefficients are defined with respect to a material axis whose orientation 
usually does not correspond to the principal axis of the structural element. The transforma- 
tion matrix used within the program to rotate the elastic coefficients was expressed ?n terms 
of direction cosines. It is very important to note that all the angles are measured from the 
local axis to the material axis (see figure 1). 


The transformation for the elastic coefficients is then 


where 


f" A 1 

— 

[VT T a 1 

"u" 

e 




L J local 

L J L J mat'l 



- 2 

2 * 
cos p 



COS CL 


r, i 

2 o 

2 


L U J = 

cos /3 

cos a 



L-2cosa cosjS 2cosa cosj3 
and for the thermal expansion coefficient vector 
a, 


cos a cosj8 i 

-cos a cos# 

2 2 P I 

cos a - cos p J 


‘11 

“22 

a 12 


= U 


-1 


V 


“22 

a 


local 


12 


mat* I 


( 12 ) 


(13) 


In the computer program all these transformations are performed automatically for each 
element. It is only necessary to indicate the direction of the material axis by specifying 
two points on the axis. Several material axes can be specified for a single structural model. 
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TWO-DIMENSIONAL ELASTIC SOLUTIONS FOR 
CRACK-TIP STRESSES AND DISPLACEMENTS 


The field equations of elasto statics appropriate to isothermal deformation under plane- 
stress conditions (a^ = T = = 0) with zero body forces are given below and represent a 

basis for all the equations developed in this section. 


e 

x 


3u 

x 

ax 


e 

/ 


au 

-r— and 

ay 



2 2 2 

a e a e ay 

X + 1 = 2Z 

iy 2 * 2 3xa y 


e 

X 


a n 

a 12 

a 16 

e 

y 


a 12 

a 22 

a 26 

Y 

xy J 


- a 16 

a 26 a 66j 



au a.u 

ay ax 



aT 

+^ = o 


ay 


and 


(14) 

05 ) 


(16) 


07 ) 


a a aT 



For plane-strain applications (e =Y = Y -0). the constants a., in (16) must be 

r rr ' Z yz ZX 'I 

replaced by their plane-strain counterparts 0.. given by 


0 .. 

'I 



a i3 q j3 

a 33 


(i/i = 1/2,6) 


08 ) 


The form of the equilibrium equations (17) implies the existence of a stress function 
U(x,y) suc h that 


CT 

X 




T 

xy 


3 2 U 


3x3y 


(19) 
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By using (19) to eliminate the stress components in (16) and substituting the resulting strain 
components into the compatibility condition (15), we find that the stress function must satisfy 

4 4 4 4 4 

cTu „ 3 U 1/n , , a u „ a u , a u_ A 

a 22 > 4 “ 2a 26 3 2 °12 a 66 * . 2. 2 ' 2a 16 3 a l 1 4 0 (20) 

ax ox ay ox oy dxoy oy 

Equation (20) may be written as the product of four first-order operators 

D rsr'*ilr = < 2 » 

If Pj, \^ 2 > ^3 anc l are roots of the characteristic polynomial, 

°] / ‘ 2 a , 6 ^ + < 2a 12 + ‘ 2 a 26 M + a 22 = 0 ® 2) 
then (20) may be written as 

D 1 D 2 D 3 D 4 U(x,y)=0 (23) 

The roots of (22) occur in complex-conjugate pairs and can be shown by energy considerations 
always to have non-zero imaginary parts. The designation u^, i-^, Pj and P 2 shall henceforth 

be used for the roots of (22), and it will be taken for granted that p. - p. 0. The general 

solution of (23) for U(x,y) can now be written by repeated integration, but the form it takes 
depends on whether or not the roots of (22) are distinct. A discussion of cases follows. 


Distinct Roots (p ^ T^Pj) 

For distinct roots, the general solution of (23) for real U(x,y) is 
U(x,y)=2Re[U 1 (z 1 ) + U 2 (z 2 )] 

where Uj and U 2 are arbitrary functions of the complex variables 
Zj = x + p^y and z 2 =x+P 2 y 

respectively. Upon substituting (24) into (19) the following expressions are found for the 
stress components. 

\ -2M^Uf(i ] )+^U“(i 2 )] 


» =2Re[UJ(z ] )+U“(z 2 )l 
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and x xy = -2Re[ l a 1 UJfzj) +^U^(z 2 )] 

in which prime (') denotes differentiation with respect to the parenthetical argument. It may 
be verified by differentiation that the displacement components corresponding to (26) are 
g iven by 

u x =2Re[p 1 U' 1 (z 1 ) +p 2 U^(z 2 )] (27) 


and u y =2Re[q 1 U' 1 (z 1 ) +q 2 U 2 (z 2 )J 


provided 


p i a n |J i +a i 2 _ a i6 |i i ' 

p 2 a n^ 2 +a 12 " °16H2 


a 22 

and q] =a ]2 ^+--a 26 , 

a 22 

q 2 a 12^2 + u 2 " a 26 

(28) 

The U.j and U 2 components of the stress function are taken to be half-power series, 

i.e.. 

® n-2 

co n-2 


u i’ (z i ) = 1 c „ z i 2 ° n,:l 

U 2< z 2> ‘ I °n z 2 2 

(29) 

n=l 

n=1 



These will now be shown to satisfy boundary conditions appropriate to free crack faces by 
properly relating the complex coefficients C and D . From (26), the stresses corresponding 
to a typical term in (29) are n n 


a 

x 


n-2 n-2 1 


= 2Re 


2 

1 





and 


a 

/ 


= 2Re 



n-2 

2 


+ D 

n 



n-2 n-2 1 


t = -2 Re 
xy 


I I-S c 
L I n 



D 

n 



(30) 


Upon considering the crack-tip neighborhood and coordinate system shown in figure 2, 
we can write 


CT J X '°) = ^ (x,0) = 0 for x <0 


(31) 
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as the crack-face boundary conditions. From (30), these require 


(C + D ) = 0 , when n is 6 X j P 

Im n n odd 


and 


Re 


(32) 


lm* 1 l C „ + ^ D n )=0 ' whenni5 odd" 


Equations (32) are satisfied if we set 


and 


C + D = (i) n+1 A 
n n n 


(33) 


In which i is the imaginary unit and and B n are arbitrary real constants. Using (33) to 

write C and D in terms of A„ and B„, we obtain 
n n n n 


,n+l ^ A n - B n 


and 


c = (i) 

n ^ “ ia l 


n+1 ~ h*"jA 

D = (i) n+1 — LH 

n 


(34) 


From (26), (27), (29), and (34), we are now able to write the crack-tip stresses and dis- 
placements corresponding to a typical term in the Uj and U 2 series. 


/«n +1 

a = 2 Re ( — ^ 

x ) \x - |ij 


n-2 n-2 x 


n-2 n-2 N 


l A n M lM*Vl "^2 Z 2 ) +B n\^2 Z 2 “^l 2 ! 


a = 2Re 

y 


jiTL 

M 2' |J 1 


n-2 n-2 N 


n-2 n-2' 


' A nh z l 2 ) +B n( z 2 2 ‘ ^ 


(35) 


and 


t = -2 Re 

xy 


= - Re 
n 


(i) n+I 




(i) n+1 


'n-2 n-2 



n-2 n-2\l 


2 2 \ / 2 2 

A n.ia. (z. - z 0 + B hUz 0 -iJ.z. 

n 1 2 1 I 2 / n\22 II 


^2-h 


Pi 1 Po 1 \ 

• A ^ r\ 


Po n p. n\ 

2 ‘ 9 \ / z 2 7 

A n (^2 2 l“ M l Z 2/ +B n( Z 2 " Z 1 


’I 


'1 


(36) 
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Repeated Roots (p.^ = P 2 = P) 

In this case the operators in (23) are not distinct; each is once repeated; i.e.. 


where 


D^U( x / y ) =° 


9 -9 , _ 9 — 9 

D =r p^— and D- p t - 

p dy 9x P dy dx 


( 37 ) 


(38) 


Repeated integration now gives 

U(x,y) =2Re[U 1 (z 1 ) + ^ U 2 (z^] (39) 

as the general real solution of (37). In (39), 

Zj =x + py and z^ =x + py (40) 

Expressions for the stress components corresponding to (39) are found from (19) to be: 

=2Re[p 2 Uy( Z] ) +2ppl)^(z ] ) +p 2 I ] U^(z 1 )] 

a y =2Re[UJ(z 1 )+2U^(z 1 )+I 1 U^'(z 1 )] (41) 

and t = -2Re[pUJ( Zl ) + (p +p) U^) +pI ] UJ(z 1 )] 

Moreover, the displacements, 

u x =2Re [pjU'jfcj) +p 2 U 2 (z ] ) +P 3 z 1 U 2 (z 1 )] 

and U y =2Re[q 1 U' 1 (z 1 ) +q 2 U 2 (z 1 ) +q 3 I 1 U 2 (z 1 )I (42) 

can be shown* to be consistent with the stresses (41), the stress-strain equations (16) and the 
strain-displacement equations (14) provided 


*ln verifying (42) subject to (14), (16), (41), and (43) certain identities must be used appro- 
priate to pand p as repeated roots of (22). These are: 

a 16 =a ll^ + ^ 2 _ 

2a 12 +a 66 = a n (p +4pp+n) 


and 


a 26 =a ll^ +|J) 

a 22 =a,/n 
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p l = P3 =li a ll +a 12“^ a 16 
p 2 =n(2fl-t4 a n +a 12 -M ]6 


22 

q 1 =q 3 = a 12^ + — ' a 26 


( 43 ) 


and q =a JI+— (2 

M 2 12 a \ \xJ 


\x/ °26 

Again, to develop a half-power series analogous to (29), we take 
“ n-2 * n-2 

u i'< z i> = 1 c „ z i 2 and U 2 (z i' = 2 D n z , 2 

n=l n=l 


(44) 


in which C n Dp are complex constants and n is anticipated to be a real integer. 
Following the same procedure used for the distinct-root case, we substitute (44) into (41) and 
consider crack-face boundary conditions (31). This leads to the requirement that C and D 
typically satisfy P 


/- J . n+ 2 ri _.n+l A 
C +—x — D = 1 A 
n 2 n n 


and uC +(V D = i n+1 B 

n \ 2 / n n 


(45) 


in which A and B are real constants. The simultaneous solution of (45) for C and D 
. , , n n n n 

yields 


.n +1 


C =- 


n U - H 


.n+1 

1 


B n] 


(46) 


and D = (B - f-iA ) 

n — n n 

u -n 

When these expressions are incorporated in (41) and (42), the following formulas for 
stresses and displacements evolve for a typical term of the series: 
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a = Re 
x 


n-2 

; n+I — 
il-U Z ' 


A U 2 L- 2\i - n(n -2)-i 
n \ z i 


+ B H (4m- - [i (n + 2) + H(n - 2) — - 

n \ 7 _ 


a = Re 


, , n-2 

i " +1 T 

? - |i 1 


A hj(n-4)+2|i-^(n-2)-l . 


1 




(47) 


and 


T = Re 
xy 


, , n-2 

• n +i — 

Z 1 

\d-\A 


^ A n< n - 2 >( ^7" y +B n (n^- 2 U-^(n -2)^1 


u = Re 
x 


i n+1 h 

Z 1 n 

l±-\l 


A n \ p ]( 2u - +n ^ " 2p 2^ “ np l^ 


+ B n( 2 P2' P l (n+2)+np li7 /J 


u = Re 

y 


i" +1 §2 

Z |n 

H -|i 


A n ( q 1 ( 2 ^ + n ^ “ 2^2^ “ nq l |a ; 


+ B n( 2q 2‘ q l (n+2)+nC, lI7 


(48) 


The repeated-root case includes the stress and displacement functions appropriate to the 
crack-tip neighborhood in an isotropic material. For an isotropic material and plane stress 
conditions. 
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II 


a ll a 22 E 


a 12 " E 


a 16 = °26 = 0 


nnd a - 2 ( ] + V ) 

and a 66 E 

With these simplifications, the characteristic polynomial (22) reduces to 
4 2 

|i +2n +1=0 

which has as a solution the repeated roots 
|i = ± i 
Consequently, 

U(x,y) = 2Re [U^ (z) + zU 2 (z)] 


in which 

z = x + iy and z = x - iy 
Upon integrating (44)^, we find 

n+2 £ 

U.(z) =C* z ^ a nd U 0 (z) = D* x' 

1 n 2 n 

as components of the isotropic form of U(x,y). Combining (52) and (54), we find 

n+2 n+2 £ ii 

U(x,y) = C* z +C* z 2 + z D* x~ + z D* x~ 

' n n n n 

In polar form (see figure 2), 


( 49 ) 


(50) 


(51) 


(52) 


(53) 


(54) 


(55) 


In (54) it is convenient to use C* and D* to represent 4C n /n(n+2) and 20^/n, respectively. 
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n n .n0 n 

2 2 ' 2 

z = r e 


n0 . . n9 
cos -j- + i sin-j 


) 


n in ._n0 n 

,-2 2 “'2 
and z = r e 


n 

__ 2 / n0 . . n0\ 

~ r V cos T " 1 s,n TV 


so that (55) becomes 


n+2 


U(r, 8) = C* r 2 (cos(j + l) 0 + i sin( £ + l) 8^ 


n+2 

— 2 
+ C* r Z 
n 


^ cos ( J + * ) 9 _ ' sin ( j + l) 9 ) 


n+2 


+ D n r 2 ('“(I -1 ) 9 +isin (f-') 8 ) 


n+2 


+ D *r 2 (cos(§- l) 9- isin(£- l) 8 


n+2 


= 2r 2 Tr«(C*) cos ( j + l) ® " lm(C*) sln( j + l) 8 
+ Re(D*)cos(j- l) 8- lm(D*)»ln(j- l) e] 


(56) 


(57) 


which is a typical term in the familiar Williams series of stress functions appropriate to the 
crack-tip neighborhood in an isotropic material (ref. 14, 15). Imposition of the crack-face 
boundary conditions (31) leads to 


+ 1 


Re(D*) = 

n n 


— Re(C*) 
n n 


y + H) 


£ +1 

2 - r'J 


( 58 ) 
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Upon substituting (58) into (57), we find the following expressions for stresses and displace- 
ments in polar form: 


a 

r 


. i au , i a u 
rar r i 2 ae 2 


:1 Re(C*) r 2 n (n +2) cos " * 2 - (n - 6) cos (^1) 8 

+ I|m(0)r 2 n^-<h+2)sln(^)e + -^2— (n - 6) sin (^) 0 


a e =' 


d 2 U 


3r 


n-2 

= i Re(C*) r 2 n 
A n 


/ _LO\ i n + 2\ n +2 
(n +2) cos ( — s — ) cos 


2 7 n + 2(-l) n 


n-2 


(59) 


n-2 

+ i lm(C*) r 2 n(n +2) 
z n 


V 2 ' n -2(-l) n 1 2 ' . 



n-2 

~ J 8®( c *) r 


n T(n +2) sin ( - * 2 ) 8 - 


+ 2 


(n - 2) sin 


n+2(-iy 


- (^) •] 


i h.K‘> ' 3 " [• <" *!) 8« (^) 8 (" ■ 8) ™ 


and 
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“r “HT Re(C ^ f2 [- <n + 2> cos (H 1 ) 9 - n " 2 ' " ~ n ( 4 -TT7 - n ) ‘“(nr) e ] 

- HT lm(C ^ f2 [ <n + 2) !in m 1 ) 9 + n "^n ( 6 - TTZ - ") !in (¥) ®J 


u e = i r iRe(c ^ r2 (" + 2 ) s ' n (^) a - 2 n ( 6 - TC7 + n ) ’ K ^) 9 

L n+2(-l) J 


_ ( 60 ) 


+ -^-=— ■ lm(C*)r 
E n 


- (n +2) cos (¥£) 0 + "J 2 (i - -p£ +n) cos (^)e 

n - 2 (— 1 ) 
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FORMULATION OF CRACKED ELEMENT STIFFNESS MATRICES 


Two types of cracked elements, as shown in figure 3, have been developed and 
implemented, because many fracture mechanics problems are symmetric about the plane of 
the crack. One formulation takes only the symmetric terms in the series solution and, 
hence, is applicable only to symmetric problems (K|| = 0); the other formulation makes use 
of both symmetric and anti-symmetric terms and is applicable to unsymmetric or mixed-mode 
problems (K| and K||). 

The coordinate system of an 8-node symmetric element has its origin at the crack tip. 

It is rectangular in shape with a three-to-one aspect ratio. Placement of the nodes relative 
to the rectangle is pre-determined with a node at each comer plus nodes at the one-third 
points of each of the long sides. The lower side (nodes 6, 7, 8, and 1) is coincident with 
the crack direction and presumed axis of symmetry. Nodes 6 and 7 are on the free crack 
face. Nodes 8 and i are on the prolongation of the crack. They are constrained rigidly as 
to vertical displacement and are free of shear forces - conditions that are consistent with 
symmetry. 

The 8-node symmetric element has 16 displacement degrees of freedom, two per node 
corresponding to the in-plane displacement components. Thus, it incorporates the first 13 
symmetric terms of the series associated with rigid-body motion in the plane. These 13 
coefficients and the three rigid-body parameters are referred to as the 16 generalized co- 
ordinates of the cracked element. The stresses and displacements corresponding to these 16 
generalized coordinates are evaluated on the boundary of the element. Products of stress 
and displacement contributing to boundary work are performed and integrated. The result is 
a homogenous quadratic form in the generalized coordinates, and the coefficient of each 
term is an element of the cracked element stiffness matrix with respect to the generalized 
coordinates. Once the stiffness matrix with respect to generalized coordinates is deter- 
mined, the stiffness matrix with respect to nodal displacements is formed using the series 
(with rigid-body terms) to write nodal displacements in terms of the generalized coordinates. 

The 10-node unsymmetric cracked element as shown in figure 3 is square with equally 
spaced nodes around its boundary. Like the symmetric cracked element, the shape of the 
element and relative location of the nodes were chosen to provide modeling convenience. 
The generalized coordinates correspond to the first 9 symmetric terms and first 8 anti-sym- 
metric terms of the series, plus the 3 rigid-body displacement parameters. The stiffness 
matrix was again generated by integration around the boundary. 

To formulate the cracked-element stiffness matrix, it is convenient to introduce 
dimensionless variables for both anisotropic and isotropic cases. 
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Anisotropic Case 


The dimensionless variables for the anisotropic case are given as 


Zj = Z ] /A , 

Z 2 = 2 2 /A 


P 1 " p l/ a ll ' 

P 2 = p 2 /a U 


Q 1 =q i/ a ll ' 

q 2 - V a l 1 


5-' 

C = a, , A A , 

nil n 

2 1 

D = a,, A B 

nil n 

(61) 

S = a. . a , 

x 1 1 X 

S *" Qi l ^ / 
y n y 

^xy a l 1 T xy 

U = u /A 
X x' 

U = u / A 

y y 


where A is a characteristic length and is taken to be the distance between nodes on the 
cracked element. In terms of the dimensionless variables in equation (61), equations (35) 
and (36) take the form below 



I 


and 


U 

x 

U 

y 


4 


n 


Re 


(O n * 1 

k 2 -h, 



Q, 


n 




+ D 

n 


Q 

2 



Qi 




( 63 ) 


Equations (62) and (63) are the basic equations needed to form the stiffness matrix of 
the anisotropic cracked element. The strain energy V stored in the cracked element can 
be computed by numerically integrating the work of the surface tractions around the bound- 
ary. This yields 



T 

q 


kq 


in which h = the uniform thickness of the cracked element 


(64) 


q = the column matrix of generalized coordinates 
q^ = the transpose of q 

k = the stiffness matrix with respect to generalized coordinates 

Let D denote the column of dimensionless nodal displacement components in the 
elemental cartesian coordinate system. The matrix C in 

D = Cq (65) 


is obtained by evaluating equations (63) at the nodes. The inverse of C may be found 
numerically with the result 


q 



D 


( 66 ) 


Thus, from (64), 


2V= D T (C' 1 ) — k C _1 D 

a ll 


( 67 ) 


or in terms of D*, the column of dimensional nodal displacements 
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2V = D* T (C -1 ) T — k C _1 D* 

°11 

L .] T 1 

K = JL. (c ') k C ' 

°11 


( 68 ) 


(69) 


Once the stiffness matrix of the cracked element is formed, its incorporation into the 
stiffness matrix of an assembly follows exactly the procedure used for conventional elements. 
The same approach was followed for the case of repeated roots. 


Isotropic Case 

The dimensionless variables introduced for an isotropic case are given below: 


R 4 ' 


s = — 

R G ' 


.* = 


2G 


U R = t 

S 

0 G 


Re (C*) 
n 




s = — 

R0 G 


(70) 


a*= - ■ 


n - 2 H) n lm(C ^ 


n 2G 

where A is also a characteristic length as defined before. 

In terms of the dimensionless variables in equation (70), equations (59) and (60) take 
the following form: 

« 5 ■ 1 

S R (R, 0) = n R | s n ( n + 2 ) cos (5 + ^ 9 + f (h) ( n " 6) COS (^ - 1) ej 

+ a* £g(n)(n + 2) sin (^ + 1)0 - (n - 6) sin (^ - 1) ©J j 
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in which 


Sg (R, 0) = ^ n (n + 2) R | s* £cos (^ + 1)0- f(n) cos - 1) sj 

+ a^* g(n) sin (^+1)0 + sin - 1) 0^ j 


n i 

00 _ — ] 

W R ' 9 > = L "R 2 {,[ 


(71) 

(n + 2) sin (^ + 1)0- f(n)(n - 2) sin - 1) 0^ 


+ a* 


g(n)(n + 2) cos + 1) 0 - (n - 2) cos - 1)0 


U r (R, 0) = 


(n + 2) cos (^ + 1) 9 - f(n)(6-8 £-n) cos (^ - 1) 9 


r %[• 

n - ! 

+ a* £(n + 2) g(n) sin (^+1)0 + (6-8 £-n) sin (^ - 1) 0 




os £ 

2 


(72) 


Ug(R, 0) = y ^ R |s* j^(n + 2) sin (^ + 1) 0- f(n) (6-8£ + n) sin (^ - 1) 0^ 
+ a* j^(n + 2) g(n) cos (^+1)9- (6-8£ + n) cos (^ - 1)0 j 


f(n) = 



j + (-!)" 


g(n) = 





for plane strain 
for plane stress 


(73) 
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Equations (71) and (72) are the basic equations needed to form the stiffness matrix of the 
isotropic cracked element. 

It is noticed that equations (72) for the displacement components have been written so 
as to keep the parts independent of £ and distinct from the parts dependent on £. This 
distinction was made to permit the storage of the stiffness matrix as the sum of two matrices 
that are each internally independent of material parameters. 

By numerical integration, the strain energy V stored in the cracked element can be 
written in terms of the stiffness matrices kj and \< 2 , i.e., 

2V = GhA 2 q T (k 1 +£k 2 )q (74) 

in which h is the uniform thickness of the element. Following the same approach as for 
the anisotropic case, equation (74) can also be written in terms of D*, the column matrix 
of dimensional nodal displacements, 

T -1 T -1 

2V = D* (C ) Gh(k ] +£k 2 )C D* (75) 


So that 


K = Gh (C 


T 

) (y^yc 


-i 


(76) 


Each of the stiffness matrices of these isotropic cracked elements is stored as the sum 
of two matrices that are independent of the size and material properties of the cracked 
element. This means that the boundary integration mentioned previously does not have to 
be repeated for each application. Consequently, the efficiency of the analysis program in 
conventional finite-element applications is not diminished at all by the addition of the 
cracked element. The displacement incompatibility that exists between nodes where the 
cracked element interfaces with a conventional element seems inconsequential in light of 
the exceptional accuracy obtained in many and varied applications, including the sample 
problems. 
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STRESS- INTENSITY FACTORS AND STRAIN-ENERGY RELEASE RATES 


Stress- Intensity Factors 

- 1 /2 

The leading terms in equations (35) and (59) contain the singularity r ; all subse- 
quent terms are non-singular. The coefficients Aj and Bj (Cj* for isotropic case) are 
related to the opening and sliding mode stress-intensity factors K| and K|| by the following 
formulas: 


Anisotropic Case : 

K | = x-” V X ' = A 1 

(77) 

K|j = x ^ m o 7 x T xy(x, o) = 2/2T B^ 

Isotropic Case : 

K, = r 'i m o /2 ttT <Jq (r, o) = 6/2T Re (C^) 

(78) 

K|l = r ^ m o v/ 2ffr T r0(r, o) = -2yJ 2n Im (Cf) 


Strain-Energy Release Rates 

The strain-energy release rates are related to the stress- intensity factors by the elastic 
coefficients in the following fashion: 

G. = cK. 2 i = I, II (79) 

where "i" indicates the mode number. The total strain-energy release rate is the summation 
of energy rate at each mode based upon linear superposition. The elastic coefficients c 
that relate energy rates to stress- in tensity factors are listed in Table I. 
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TABLE 


ELASTIC COEFFICIENTS RELATING ENERGY RATES TO STRESS-INTENSITY FACTORS 

(Reference 16) 


MATERIAL 

CONDITION 

Isotropic 

Plane-Strain 

Plane-Stress 

Ortho tropic 

Plane-Strain 

Plane-Stress 

Anisotropic 

Plane-Strain 

Plane- Stress 


ELASTIC COEFFICIENTS 


MODE I 



°H a 22 



Im 

("1 + »*2)] 


^ Im 

[ a ll 1 + **2> 





































BRIEF DESCRIPTION OF CRACKED ELEMENT COMPUTER PROGRAM 


The computer program is an a 1 1- FORTRAN program which uses only conventional l/O 
(FORTRAN READ and WRITE statements) and the standard library routines. It contains the 
following finite elements: 

• Axial element 

• Linear spring element 

• Plane-stress/ strain, isotropic and anisotropic, triangular element 

• Eight-node symmetric isotropic cracked element 

• Eight-node symmetric anisotropic cracked element 

• Ten-node unsymmetric isotropic cracked element 

• Ten-node unsymmetric anisotropic cracked element 

The computer program consists of a main program called CRAK, which allocates the 
core and calls sequentially four main driver subroutines (INPUT, KFORM, CHOL, 
OUTPUT), and 32 subroutines. These subroutines have been overlayed to reduce the core 
required for program code, permanent data, and program table. The segment structure of 
the program is shown in figure 4. The name and function of each subroutine are given in 
Table II. 

The input segment of the program contains data generator features to minimize the 
labor required to input a model. During the input sequence, a banding algorithm is used 
to resequence the node numbers to minimize the band-width of the structural stiffness 
matrix which is assembled in the KFORM segment. This banding algorithm is quite effec- 
tive and helps to minimize execution times. A banded Cholesky decomposition procedure 
(segment CHOL) is used to solve the system of equations. The computer program organiza- 
tion is especially characterized by its flexibility of application. The planned segmented 
structure of the program also makes it easy to insert additional cracked and conventional 
elements as desired. 

The scheme of the analysis is basically divided into the following stages: 

• Input of geometry, material data and restraints 

• Input of applied nodal loads and imposed nodal displacements 

• Calculation of element stiffness matrices 

• Assembly of element stiffness into a system stiffness 

• Solution of simultaneous equations for nodal displacements 

• Output nodal displacements 
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• Calculation and output of forces (or stresses) in conventional elements 

• Calculation and output of stress-intensity factors and strain-energy release 

rates 

The overall logic flow of the computer program is depicted in figure 5. 

The computer program is efficient in its use of computing time and core storage. Core 
use can easily be changed to fit individual problem sizes if desired. Full advantage has 
been taken of symmetry and the banded nature of the simultaneous equations to be solved, 
and no data are generated when not required. 
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TABLE II 


FUNCTION OF SUBROUTINES 


SUBROUTINE 

NAME 

FUNCTION 

INPUT 

Inputs model general specifications 

CORIN 

Inputs nadal coordinates and restraints 

TRIN 

Inputs triangular elements for both isotropic and anisotropic cases 

SPRGIN 

Inputs spring elements 

AXIN 

Inputs axial elements 

CRKIN8 

Inputs 8-node symmetric cracked elements for both isotropic and 
anisotropic cases 

CRKN10 

Inputs 10-node unsymmetric cracked elements for both isotropic and 
anisotropic cases. 

MATIN 

Inputs material properties 

DISPIN 

Inputs imposed displacements 

LOAD IN 

Inputs applied loads 

CONNEC 

Generates data which subroutine BANDIT requires 

BANDIT 

Renumbers the node points to minimize band width 

BAND 

Determines the band-width 

KFORM 

Forms a complete stiffness matrix 

TRI 

Computes triangular element stiffness matrix 

SPRNG 

Computes spring element stiffness matrix 

AXI 

Computes axial element stiffness matrix 

CRAK8 

Switches control to either CRAK8I or CRAK8A 
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TABLE II (Continued) 
FUNCTION OF SUBROUTINES 


SUBROUTINE 

NAME 

FUNCTION 

CRAK8I 

Computes stiffness matrix for 8-node symmetric isotropic cracked 
element 

CRAK8A 

Computes stiffness matrix for 8-node symmetric anisotropic cracked 
element 

CRAK10 

Switches control to either CRK10I or CRK10A 

CRK10I 

Computes stiffness matrix for 10-node unsymmetric isotropic cracked 
element 

CRK10A 

Computes stiffness matrix for 10-node unsymmetric anisotropic cracked 
element 

INVERT 

Inverts a matrix 

CPLXRT 

Computes complex roots of a specified 4th order characteristic polynomial 

CPWR 

Raises complex number to a real power 

THCOEF 

Transforms thermal coefficients from one reference axis to another 

THDISP 

Computes thermal displacements for both symmetric and unsymmetric 
cracked element 

CHOL 

Solves simultaneous equations 

OUTPUT 

Controls output data 

DISPOT 

Outputs Displacements 

TRIOUT 

Calculates triangular elements output 

SPROUT 

Calculates spring elements output 

AXIOUT 

Calculates axial elements output 
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TABLE II (Continued) 
FUNCTION OF SUBROUTINES 


SUBROUTINE 

NAME 

FUNCTION 

CRK08 

CRKOIO 

Calculates Kj for the 8-node symmetric elements 
Calculates Kj and Kjj for the 10-node unsymmetric elements 


32 









SAMPLE PROBLEMS 


This finite-element computer program has been tested extensively. Both symmetric and 
unsymmetric cracked elements have performed well with respect to accuracy and efficiency 
as shown by sample problems presented here. These results were achieved with a relatively 
coarse finite-element grid. With refinements in the grid, even more accurate results would 
be obtained. The symmetric cracked element normally gives more accurate results than the 
unsymmetric cracked element, which is understandable in view of the fact that the unsym- 
metric cracked element has fewer degrees of freedom that it can bring to bear on the first 
mode. However, the unsymmetric cracked element can be used in a much wider class of 
crack problems and is more practical for industrial applications. 

To illustrate the capabilities of these two types of crack elements (symmetric and un- 
symmetric), a number of selected sample problems are included in this report. These sample 
problems were chosen to demonstrate (1) the accuracy and economy of the elements, and 
(2) the versatilities of the elements to perform analyses for structural configurations of prac- 
tical importance. 


Anisotropic Case 


Center- Cracked and Double-Edge-Cracked Qrthotropic Tension Plates - A center- 
cracked and/ or a double-edge-cracked orthotropic plate (Figure 6) subjected to uniform 
tension as solved by Snyder and Cruse (ref. 17) was investigated. The plate analyzed is 
228.6 mm long and 76.2 mm wide. The finite-element models used for both the center- 
cracked plate and the double-edge-cracked plate are shown in figure 7. The half-crack 
lengths, a, analyzed are 15.24 mm and 22.86 mm for the center- cracked plate and 
15.24 mm only for the double-edge-cracked plate. 


Numerical results were obtained for 10 representative graphite fiber-reinforced epoxy 

laminates: (0) , (±30) (±34) , (±45), (±56), (±60) (90) (90 /±45) , (90 /±45) , 

/ . s . s s. s s s s 4 s z s 

(0/—45/90) . Lamina properties used were: 


E n = 144.795 GPa , 


”* 9.653 GPa 


E 22 = 11.722 GPa 


= 0.21 


Table III indicates the complex roots obtained from the characteristic equation for each 
laminate. 

Results are presented in Tables IV, V, and VI. The distribution of differences from 
Snyder and Cruse's results are summarized in table VII. Notice that the results generated 
using 8-node symmetric cracked element are closer to Snyder and Cruse's results than these 
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TABLE III 


COMPLEX ROOTS OF CHARACTERISTIC EQUATION 


MATERIAL 


^2 

“l 

* 

a 2 

■ 8 2 

(0/i 45/90) s 

0 . 

0.9999 

0 . 

1 .0001 

(90 2 /t 45) $ 

0.2680 

0.7007 

-0.2680 

0.7007 

(90 4 ^ 45) s 

0 . 

0.5554 

0 . 

0.8372 

(90) $ 

0 . 

0.2704 

0 . 

1.0522 

(±45) s 

0.7763 

0.6304 

-0.7763 

0.6304 

(- + 30) $ 

0.9648 

1 . 0090 

-0.9648 

1 .0090 

(i60) s 

0.4951 

0.5177 

-0.4951 

0.5177 

<°>S 

0 . 

0.9504 

0 . 

3.6982 

< ±34 >S 

0.9415 

0.8730 

-0.9415 

0.8730 

(±56) s 

0.5711 

0.5296 

-0.5711 

0.5296 


NOTE: /i. - a. + j8. I , where i - V-T" , j - 1 , 2 


























































TABLE IV 


STRESS-INTENSITY FACTORS OF 
CENTER-CRACKED PLATES (a = 15.24 mm, a/w =0.4) 


MATERIAL 

SNYDER & CRUSE 
K j (MPa - /m) 

8-NODE 

CRACKED ELEMENT 

10-NODE 

CRACKED ELEMENT 

K . 

% OFF 

K i 

% OFF 

(0/- 45/90) ^ 

1.670 

1.675 

40.32 

1.622 

-2.88 

PO/t 45 ) s 

1.679 

1.696 

+1.01 

1.652 

-1.61 

P0/i 45) s 

1.670 

1.702 

+1.92 

1.671 

+0.06 

P0) s 

1.683 

1.669 

-0.83 

1.755 

44.28 

^45) s 

1.745 

1.634 

-6.36 

1.617 

-7.34 

(' 30) S 

1.712 

1.643 

-4.03 

1.574 

-8.06 

(“ 60) s 

1.719 

1.683 

-2.09 

1.649 

-4.07 

(°) s 

1.645 

1.652 

+0.43 

1 .578 

-4.07 
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TABLE V 


STRESS-INTENSITY FACTORS OF 
CENTER-CRACKED PLATES (a =22.86 mm, a/w =0.6) 


MATERIAL 

SNYDER & CRUSE 

8-NODE 

CRACKED ELEMENT 

10-NODE 

CRACKED ELEMENT 


Kj (MPa - y/m ) 

K i 

% OFF 

K i 

% OFF 

(0/- 45/90) s 

2.396 

2.396 

40.01 

2.326 

-2.92 

(90^i 45) s 

2.422 

2.444 

40.91 

2.389 

-1.36 

(90^/- 45) $ 

2.396 

2.446 

4-2.09 

2.413 

40.71 

(9°) s 

2.394 

2.419 

+1.04 

2.574 

+7.52 

(- + 45) s 

2.599 

2.421 

-6.85 

2.393 

-7.93 

i 

(- + 30) s 

2.518 

2.405 

-4.49 

2.302 

-8.58 

60) S 

2.533 

2.473 

-2.37 

2.430 

-4.07 

(0) s 

2.314 

2.336 

40.95 

2.229 

-3.67 
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TABLE VI 


STRESS- INTENSITY FACTORS OF 
DOUBLE-EDGE-CRACKED PLATES (a = 15.24 mm, a/w = 0.4) 


MATERIAL 

SNYDER & CRUSE 

8-NODE 

CRACKED ELEMENT 

10-h 

CRACKED 

'JODE 

ELEMENT 


K, (MPa - V m ) 

K i 

% OFF 

K , 

% OFF 

(0 /- 45/90) s 

1.722 

1 .701 

-1.22 

1.649 

-4.24 

(90^ 45) s 

1.722 

1.729 

-HD. 41 

1.685 

-2.15 

< 50 /- +45 >s 

1 .706 

1.725 

+ 1.11 

1 .638 

- 1.06 

(90) s 

1.709 

1.689 

-1.17 

1.752 

+2.52 

( i45 )s 

1.700 

1.757 

+3.35 

1.766 

+3.88 

(-30) s 

1.783. 

1.722 

-3.42 

1.697 

-4.82 

(“ 60 ) s 

1.717 

1.767 

+2.91 

1.736 

+1.11 

<°>s 

! 1.670 

i 

1.660 

-0.60 

1.598 

-4.31 

(-34) s 

1.750 

1.737 

-0.74 

1.724 

-1.49 

(-56) s 

1.745 

1.765 

+1.15 

1.746 

+0.06 
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TABLE VII 

DIFFERENCES FROM SNYDER AND CRUSE'S RESULTS 




DISTRIBUTION OF DISCREPANCY 

TYPE OF PLATE 

TOTAL NUMBER 
OF PROBLEMS 

DIFFERENCE 

NO. OF PROBLEMS 



p/o) 

8-NODE CRACKED EL. 

10-NODE CRACKED EL. 



0 -±3 

6 

3 

CENTER-CRACKED 

8 

±3 - ±5 

1 

3 

(a = 15.24mm, a/w = 0.4) 


±5 - ±7 

1 

0 



±7 - ±9 

0 

2 



0 - ±3 

6 

3 

CENTER-CRACKED 

8 

±3 - ±5 

1 

2 

(a = 22.86 mm, a/w * 0.6) 


±5 -±7 

1 

0 



. ±7 - ±9 

0 

3 



0 - ±3 

8 

6 

DOUBLE-EDGE-CRACKED 

10 

±3 - ±5 

2 

4 

(a * 15.24 mm, a/w = 0.4) 

i 

±5 - ±7 

0 

0 

j 


±7 - ±9 i 

0 

1 

0 























obtained from the 10-node unsymmetric cracked element. Since Snyder and Cruse used the 
boundary- Integra I approach, no definite conclusion can be drawn with regard to the abso- 
lute accuracy of their results. Besides, these finite-element solutions used a large integra- 
tion step size (10 steps between nodes) along the boundaries for the numerical integration in 
forming the stiffness matrix of the cracked element. Accuracy would be improved by re- 
ducing the integration step size (increasing the number of steps between nodes). 

Longitudinally Cracked Orthotropic Strip - The longitudinally cracked orthotropic strip 
problem as shown in figure 8 was solved under plane-stress and imposed displacement condi- 
tions. An analytical solution was developed based upon the energy equivalence and Sih's 
results of ref. 18. The solution is given below as 



where a j | = Elements of the material compliance matrix 

5 = Imposed displacement 

h = Height of the plate 

Numerical results were obtained for the two different laminates shown in table VIII. 

In one case, the laminate is stiffer in the direction normal to the crack. In the other case, 
the laminate is-stiffer in the direction of the crack. Both laminate cases were solved using 
a single cracked element and the two different model representations shown in figure 9. 

The results are summarized in table IX. 

45-Degree Cracked Finite Orthotropic Plate - This problem is shown in figure 10. The 
finite-element model (figure 11) representing this plate consists of 130 nodes, 210 aniso- 
tropic triangular elements, and two 10-node anisotropic cracked elements. The problem 
was solved under uniform remote stress and plane-stress conditions. The stress- intensity 
factors obtained are given below: 

K | = 1.1 376 x 1 0 ^ MPa - y m and ^ 1 1 = ^ ^ x ^ ^ MPa - y m 

Both K| and K|| are within percent (-1 .77% and +1 .75%, respectively) of Sih's solution, 
ref. 18, (K| and K|| = 1.1581 x 10“^ MPa - fm) for an infinite sheet, and they satisfy 
NASA's requirement for the accuracy of the program. 
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TABLE VIII 


LAMINATE STIFFNESS MATRICES (Scotchply 1002) 


LAMINATE 1 

LAMINATE II 

[A] = 
(GPa) 

”ll.500 0.572 0 

0.572 34.503 0 

_ 0 0 4.854_ 


[A] = 
(GPa) 

34.503 0.572 0 

0.572 11.500 0 

0 0 4.854_ 



NOTE: 



■ M l«l*y 


where x-y refers to the reference coordinates 




TABLE IX 


SUMMARY OF RESULTS FOR A LONGITUDINALLY 
CRACKED ORTHOTROPIC STRIP 


LAMINATE I 


Analytical Solution 
Ref. (18) 




LAMINATE II 


(MPa- \J~m ) 


12.477 


ERROR 

% 


ERROR 


(MPa- v/TrT) 



1 Cracked Element 


13.845 10.96 5.515 0.75 


1 Cracked Element & 
8 Triangular Elements 


1 Cracked Element & 

1 17 Triangular Elements 


13.527 


12.474 



0.97 


-2.70 























Isotropic Case 


Cracked Tension Plates - The sample problem, shown in figure 12, was the first- 
problem analyzed with the symmetric cracked element. It exhibits the degree of accuracy 
which has been consistently achieved in numerous subsequent problems. The finite-element 
model has 30 nodes, 35 triangular elements, and one 8-node symmetric cracked element. 
The three configurations (the single-edge crack, the double-edge crack, and the center- 
crack) were all individually analyzed with this one model for an a/w ratio of 1/3. The 
model grid, which is quite coarse, results in the single - edge-crack model having 57 dis- 
placement degrees-of-freedom (DOF), while the double-edge-crack and center-crack 
models have 51 DOF. The stress-intensity factors computed using these configurations are 
compared with ASTM (American Society for Testing and Materials) values. The accuracy 
of the finite-element predictions are impressive (<1 .5% error) for all three cases. Refine- 
ments in the model would produce steady convergence toward ASTM values. 

The same cracked problems were solved using an unsymmetric cracked element. The 
finite-element model as shown in figure 13 consists of 54 nodes, 69 triangular elements, 
and one 10-node unsymmetric cracked element. The results are not as good as those ob- 
tained using the 8-node symmetric cracked element, because the 8-node element has more 
degrees of freedom to bear on the first node. 

An eccentric crack problem as shown in figure 14 was also solved, and the results 
compared with Isida's solution. 

All these results are summarized in table X for reference. 

Bi-Material Cracked Plate - This problem, as shown in figure 15, was solved in 
response to the NASA's Request for Proposal (RFP 1-31-4957), dated 28 May 1974. It was 
required that the stress-intensity factors of this isotropic, generalized plane-strain problem 
must be within i2% of Erdogan's results for an infinite plate (ref. 19). All dimensions and 
material properties are given in figure 15. 

In the finite-element model shown in figure 16, only the upper half of the plate has 
been modeled due to its symmetry with respect to the x-axis. It contains 100 nodes, 148 
plane-strain triangles, and 2 eight-node symmetric cracked elements. 

A uniform displacement of 46.228 x 10 ^mm was imposed on the top boundary nodes. 
This was arrived at as follows. For plane strain, 

H(1 -v, 2 ) S, H(1 -v 2 ) S 9 

6 = 1 = “ (81) 

C 1 2 

where H is the half-height of the finite sheet. was taken as unity. 
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TABLE X 


SUMMARY OF TYPICAL RESULTS FOR ISOTROPIC TENSION PANELS 


CRACK 

DESCRIPTION 

EIGHT-NODE 

SYMMETRIC ELEMENT (ESE) 

TEN-NODE 

UNSYMMETRIC ELEMENT (TUE) 

F. E. MODEL 

K % ERROR 

F. E. MODEL 

K ( % ERROR 

SINGLE EDGE CRACK (a/W = 1/3) 

35 TRIANGULAR 
ELEMENTS PLUS 

l 

-1.2 

1 

) 

68 TRIANGULAR 
ELEMENTS PLUS 

-2.6 

DOUBLE EDGE CRACK (a/W = 1/3) 

1 ESE, 

-1.5 

1 TUE, 

-2.6 

CENTER CRACK (a/W = 1/3) 

DOF =57 

+0.15 

DOF = 105 

-1.2 

ECCENTRIC CRACK 
(IS IDA'S PROBLEM) 

116 TRIANGULAR 
ELEMENTS PLUS 

-1.1 

(SIDE A) 

NONE 

NONE 


2 ESE, 

+0.5 




DOF = 152 

(SIDE B) 



















The computed stress- in tensity factors for the two crack tips are 

K ] = 2.67496 x l(f 3 MPa -/nT K. 2 = 0.25067 x 10" 3 MPa -/i iT 

These values include in their definition a factor of yfrr~. To compare with the values in 
the Request for Proposal, it is necessary to divide by^/lr"; hence, 

K ] = 1. 50918 x 10~ 3 MPa -/nT K 2 = 0. 14143 x 10 -3 MPa -/nT 

These values are quite close to those computed by Erdogan and Biricikoglu in ref. 19. The 
percentage of error is -0.12% and ■HD. 75%, respectively. The axial stresses and displace- 
ments along the crack line are plotted in figures 17 and 18. 

Stiffened Plate - The rate of growth of a crack in a stiffened plate will be influenced 
by the presence of the risers. The simple stiffened plate configuration shown in figure 19 
has been analyzed to determine this effect. The plate is loaded by a stress along its edge 
normal to the crack. Two cracks are assumed to propagate from the center of the two edges 
of the plate toward the risers. When the cracks have passed underneath the risers, they are 
assumed to extend up the riser and into the plate at an equal rate until the riser is com- 
pletely failed. The crack is then assumed to continue to propagate into the plate. 

Only one quarter of the configuration shown in figure 19 is idealized, the remainder 
being represented by symmetric boundary restraints. Both the plate and the riser are 
idealized using triangular membrane elements. The 8-node crack element is used to repre- 
sent the crack tip. 

The results of the analysis is shown in figure 20 in the form of the crack tip symmetric 
stress- intensity factor plotted against the crack length. It is evident that the riser will 
have a considerable effect on the rate of crack growth. 

Thermal Stress- Intensity Problem - Sih (ref. 21) has given the following plane-stress 
formulae for the stress-intensity factors appropriate to a crack with constant-temperature 
faces (see figure 21) acting as sinks for the steady, radially inward flow of heat Q at 
infinity. 

EorQO+l')/^- and K n = O (82) 

1 4^k " 

In (82), a is the coefficient of thermal expansion, and k is the thermal conductivity. 

We now solve the associated temperature-distribution problem to obtain thermal input in a 
form acceptable to the computer program. 

The general real solution for the temperature, T, in a steady heat-conduction problem 
in two dimensions is 
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T(x, y) = f(z) + f(z) 


(83) 


where f Is an analytic function of z = x + iy. The symmetry of the problem indicated in 
figure 21 permits us to write 


T(x, y) = T(x, -y) 
which leads at once to 

f(z) = f(z) 

permitting (83) to be written as 
T(x, y) = f(z) + f(z) 

The constant-temperature condition for the crack faces requires that 

g = f'(z) + f© 

vanish on the crack faces; i.e. 

f' + (x) + f 1 (x) = O , x 2 < a 2 


(84) 


(85) 


( 86 ) 


(87) 


( 88 ) 


In (87) and (88), prime (') denotes differentiation with respect to the parenthetical variable, 
while the superscript (+) or (-) indicates a value taken in the upper and lower half planes, 
respectively. 


The boundary condition at infinity is 

•2n 


<J ar 

* e\ 


rd 0 - Q 


(89) 


in which r and 0 (shown in figure 21) are defined by 
i0 


z = re 


(90) 


-1 


In view of (89), f'(z) is expected to vary at infinity like z 

We now make the usual substitution for a cut-condition like (88); i.e., 
g(z) =\/z 2 - a 2 f'(z) , 


(91) 
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z - a is that branch of (z - a ) varying like z at infinity. Because of 

the anticipated remote behavior of f'(z), we expect g(z) to tend to a constant value at 

infinity. 


The problem case in terms of g(z) is then 

g + (x) - g (x) = O , x 2 < a 2 ; (92) 

x “". ) ■ «, < 93 > 

The obvious solution being 

g(z) = c 1 (94) 

which leads to 

f(z) = c 1 in (z +^z 2 - a 2 ) + c 2 (95) 


when (91) is integrated. 


The constant C 2 in (95) only sets the temperature reference, which is thus far unspec- 
ified. Taking - 0 with no loss in general applicability, we find 


ST _ 

dr r 


(96) 


for large z. Upon imposing (89), we find 

= Q 

C 1 4kff 

Thus, from (86), (95), and (97), 

T(x, y) = |^in (z + ^/z 2 - a 2 ) + in (z +^z 2 - a 2 )J 

and the constant crack-face temperature is given by 

T = Q* n a 
c 2kir 


(97) 


(98) 


( 99 ) 


Referring to figure 22, it is convenient to let 
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r j = V(x + a) 


2 . 2 

+ y 


0. = tan ^ ^ — 

1 x + a 


r x 


2 . 2 

+ y 


= fnn ' 1 I 


0 — tan 


= V( X " a > 


2 . 2 

+ y 


-l 


0 O = tan 

2 x - a 


then the temperature distribution can be written as 


0 , + 0 « 2 


T(x, y) = in | (r cos 0 +/rj r 2 cos - 2 2 ) 


0 , + 0 2 2 

+ (r sin 0 + x /r^ ^ sin — ) 


(100) 


Figure 2 3 shows a finite-element representation of the first quadrant of the problem in 
figure 2 1 . The eight-node symmetric cracked element ABCD represents the crack-tip 
neighborhood in the finite-element model. Input particulars are taken to be 

Y = 255.928°K and a=114.3mm (101) 

leading to a crack-face temperature from (99) of 

T c = 255.505°K (102) 

and triangular-element temperatures (centroidal) from (98). 

The computer program was executed twice with additional materfal properties: 

E = 68.95 GPa , t/ = 0 . 3 and a = 1.8x 10 ^ m.m ^.K (1 03) 

In the first execution, the temperature of the cracked element was taken to correspond 
uniformly to that of the crack faces (102). The computed stress-intensity factor 

Kj = 4.9137 x 10" 3 MPa (104) 

is 14.95% higher than the analytical value of 4.2747 x 10 MPa This is not 

unexpected in view of the fact that the entire cracked element is taking the crack-face 
temperature, a condition thermally more severe than the analytical temperature distribu- 
tion (98). This is reflected in the results of the second execution, in which the tempera- 
ture of the cracked element was taken to be 255.567°K, which corresponds to the average 
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of six superimposed triangles (shown dashed in figure 23). The stress- intensity factor 
computed for the second execution was 

K ( = 3.6325 x 10“ 3 MPa -/7T (105) 

which is 15% lower than the analytical value. It is clear that more refinement in the 
crack-tip neighborhood would lead to a more tolerable discrepancy. It is important to 
point out that such refinement will not be necessary in more routine applications where the 
temperature gradients near the crack tip are less severe. 
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APPENDIX - ELEMENT STIFFNESS MATRICES 


Axial Element 


The positive sign convention for the axial element is illustrated in figure A-l . 
this convention k and P ^ are as follows: 




- AEaT 



For 


Triangular Membrane Element 

The sign convention for the triangular membrane element is illustrated in figure A-2. 
Note that six displacements or degrees of freedom are present. It is then reasonable to 
assume the following displacement field. 

u = a^ + 02 ^ + a^y v = a^ + a^x + a^y 

The displacements of the element can then be expressed in terms of the a. coefficients as 
follows: 


' U 1 


1 


h\ 

u 2 


1 X 2 


a 2 

u 3 


CO 

X 

CO 

X 


a 3 

< 

V 1 


1 

1 

04 

v 2 


1 x 2 


a 5 

\ V 3J 


I 

X 

CO 

X 

CO 

1 


a A 

i 6 / 


or in matrix notation 

{•I ■ H M 

then 

W ■ M ' (•) 
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where 


!>]*' ■ 


x 2 y 3 


x 2 y 3 

" y 3 y 3 


x 32 " x 3 x 2 


x 2 y 3 

-y o 


32 


where x ^ ~ ' x 2‘ 


The strains can also be expressed in terms of the displacement coefficients: 


€ = 
x 


5u 

bx 



3v 

y oy 6 


b u ov 

y - — + — 

xy oy cx 


a 3 + a 5 


Then 


\ 


— -1 

n } n n n 0 


fo } \ 

€ 

X 




a 2 

j € 

y y j 

> m 

0 0 0 0 0 1 

. 

1 a 3\ 
) a4 ( 

y 


0 0 10 10 


a 5 

\ xy i 


— — 


\ 6 / 

M-HW 



If the strains and stresses are constant within the membrane element, the strain energy can 
be expressed as 



where A Q is the area of triangular element. Let the stress-strain relations for the anisotropic 
material in the elemental reference system be written as 
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W-Mtt 

Then the strain energy can be expressed in terms of the node point displacements as follows: 

-JvM'WHMH 

And the j^kjmatrix for the triangular membrane element is, therefore, 

TnT 


|Y|= A Q t C _1 B A B C " 1 


If the multiplication is carried out and if the elements of the matrix are arranged so that the 
convention in figure A-3 applies then the kj - elements of the stiffness matrix are 

(—2j k 11 = y 3 A 11 “ 2y 3 x 32 A 13 + *32 A 33 

(-^°j k \2 = y 3 A 13“ y 3 x 32 ^ A 12 + A 33^ +X 32 A 23 


4 ^) k l 3 = ' y 3 2 A 1 1 + y 3 (x 3 + x 32 ) A 1 3 “ x 3 x 32 A 33 


(^°) k 14 = " y 3 2 A 13 + X 3 y 3 A 12 + y 3 x 32 A 33 " X 3 x 32 A 


23 



(“•) 


51 k 15 x 2 y 3 A 13 + x 2 x 32 A 33 


k 16 “ x 2 y 3 A 12 + x 2 x 32 A 23 


k 22 y 3 A 33 " 2 y 3 x 32 A 23 + x 32 A 22 


k 23 = " ^ A 1 3 + x 3 y 3 A 33 + y 3 x 32 A 12 ~ x 3 x 32 A 23 
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" y 3 A 33 + y 3^ x 3 + x 32^ A 23“ X 3 X 32 A 22 

= x 2 x 32 A 23 ‘ x 2 y 3 A 33 
= x 2 x 32 a 22 - * 2 y 3 a 23 

= y 3 2 A 1 1 ’ 2x 3 y 3 A 13 +x 3 2 A 33 

= y 3 A 13 ’ x 3 y 3 ^ A 12 + A 33^ +x 3 A 23 

= x 2 y 3 A 13 -x 2 x 3 a 33 

= x 2 y 3 A 12 - x 2 x 3 A 23 

= y 3 A 33" 2x 3 y 3 A 23 + x 3 A 22 
= -x 2 x 3 A 23 + x 2 y 3 A 33 

= -x 2 x 3 A 22 +x 2 y 3 A 23 

= x 2 A 
2 33 

= x 2 A 23 

= x 2 A 
2 22 
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where A q = ^ x 2 y 3 


x 32 x 3 " x 2 


The values of P r( . for the triangular plate element can be computed in terms of the 
element stiffness matrix. Note that 

M -MM * W 


If { P } = 0 and { p } is selected as the thermal displacement for free expansion, then the forces 
for complete restraint are 

M-MU 

For the triangular plate element, the components of jp \are 



— n/* 


< a *l 1 T x 2 

Tx,+0i* 2 Ty 


Ml 

v,=0 
v 2 = 0 

v 3 = a 22 T y 3 


The symbol * implies that the thermal expansion coefficients have been rotated into the 
local reference axis for the element. 


Spring or Fastener Element 

Figure A-4 illustrates two plates connected by a shear pin, and the direction vectors 
for the system of forces that are assumed to act on the pin are illustrated in figure A-5. For 
the spring element, it is necessary to define four node points. Node points (1) and (2) define 
the element itself. In most structural models these node points may be coincident before any 
strains occur, that is, (Z] - Z 2 ) = 0. Node points (3) and (4) are used in conjunction with 
node point (1) to define direction vectors for spring elements K] and l< 2 , respectively. K 3 
is assumed to be a spring element perpendicular to the plane containing Kj and K^. In 
other words, the system of spring elements is that illustrated in figure A- 6 . If we let u, v, 
w, be displacements, the element forces for the spring system are 
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P xl = K 1 (u l ' u 2> 


P yl “ K 2 <»1 - *2> 
P zl = K 3 (w l ' w 2> 
P x2 = ~ P xl 
P y2 = ‘ P yl 

P z2 = ' P zl 


Hence, the stiffness matrix is 


K ! 

0 

0 

- K 1 

0 

0 


0 

K, 

A 

0 

0 

-K„ 


0 

0 

K, 

V 

0 

0 

-K, 


-K ] 0 0 

o -k 2 0 

0 0 - K 3 

K 1 0 0 

0 K 2 0 

0 0 K s 


where Kj, K 2 , and are linear spring stiffnesses for the fastener element. 


10- Node Cracked Element Thermal Effects 

Let the coordinate system for the ten node cracked element be the one illustrated in 
figure A -7. The stress-strain law can be written as 


/ € 


I X 


J e 


I y 


y . 


V x/i 



a ll a 12 a 16 


a 12 °22 a 26 



a \ 

fa, ,T \ 

1 

x 1 1 

1 1 1 


M + 

Q! 22 T j 

1 

. r J l 

a 12 T ) 
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Assume a constant temperature distribution over the element; then, for free expansion 


crV = 


and 

_ 5u _ 

€ = r- = a. ,T 
x ax 11 

2>v _ 

V = 57 “22 T 

r = P- * P- = a „T 

xy ay ax 12 

The displacement function is then of the form 
u =a n Tx + C 2 a ]2 Ty + Cj 

v = a 22 Ty + C 4 a 12 Tx + C 3 

If u^ = v.| = = 0, the displacement functions become 


u =Q' 1 ^x +a 12 Ty 


v = a 22 Ty 


As in the case of the triangular plate element 


rt 


where the element of p are 

r o 


o 

II 

< 

1 

II 

CO 

> 

O 

II 

“4 = ' a 12 T4 

u 2 = -a^ T A 

< 

i 

ii 
> 

v 2 = 0 

u 5 = a ll T4 "“ 

u 3 = '“n Ti ‘“i2 TA 

v 5 = -“22 TA 


12’ 


55 


I 



u 6 = “ll Ti 

v 8 = “22 TA 

o 

II 

NO 

> 

u 9 = _<i n TA +a i2 T A 

u 7 - jTA + ^ 12 T A 

v 9 =a 22 T4 

v 7'“22 T4 

c 

o 

II 

J5 

— H 

I> 

u 8 =a 12 TA 

< 

O 

II 

o 
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til III 

! LEVEL o! LEVEL 1 | LEVEL 2 [ LEVEL 3 1 LEVEL' 4 ! 


CRAK 

INPUT 


KFORM 

CRAK8, CRAK10 


CORIN, DISPIN, LOADIN, MATIN 


TRIN, SPRGIN, AXIN, CRKIN8, 
CRKN10, CONNEC 


BAND, BANDIT 


TRI 


SPRNG 


AX I 

CRAK8I 


CRAK8A 


CRK10I 


CRK10A 


CHOL 

OUTPUT 


NOTE: 

(I) - ISOTROPIC 
(A) - ANISOTROPIC 


PI SPOT 
TRIQUT 

SPROUT 

AXIOUT 

CR KOT3 

CRKOIO 


THCOEF (A) 


THCOEF (A) 
INVERT (I, A) 


THDISP (I, A) 


CPLXRT (A) 
CPWR (A) 


THCOEF (A) 


FIGURE 4. SEGMENT STRUCTURE FOR COMPUTER PROGRAM 



FIGURE 5. OVERALL LOGIC FLOW OF COMPUTER PROGRAM 













CENTER CRACK DOUBLE-EDGE CRACK 



OOQ 

6 mm 


— | 30.48 

JLA o . 


h i ■ 

(or 45.72 mm) 

1 

r 

15.24 mm 15.24 mm 










FIGURE 6. CENTER-CRACKED AND DOUBLE-EDGE-CRACKED 
ORTHOTROPIC TENSION PLATE 
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8- NODE SYMMETRIC 
CRACKED ELEMENT 



hC RAC H 


10-NODE UNSYMMETRIC 
CRACKED ELEMENT 



FIGURE 7. CRACKED ORTHOTROPIC TENSION PLATE 
FINITE-ELEMENT MODELS 
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FIGURE 9. FINITE-ELEMENT MODELS FOR A LONGITUDINALLY CRACKED ORTHOTROPIC STRIP 









FIGURE 10. 45° CRACKED FINITE ORTHOTROPIC PLATE 





FIGURE 11. 45° CRACKED ORTHOTROPIC PLATE 

FINITE-ELEMENT MODEL 
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Kj % ERROR = -1.2 





FINITE ELEMENT MODEL = 57 DOF, 35 CONSTANT-STRAIN TRIANGLE 



8-NODE SYMMETRIC ELEMENT 


FIGURE 12. 

o 

o 


SYMMETRIC CRACKED ELEMENT MODEL FOR 
SINGLE-EDGE, DOUBLE-EDGE, AND CENTER 
CRACKED TENSION PANELS (a/w = 1/3) 




DOF = 105 


ELEMENTS = 69 


FIGURE 13. UNSYMMETRIC CRACKED ELEMENT MODEL FOR 
SINGLE-EDGE, DOUBLE-EDGE, AND CENTER 
CRACKED TENSION PANELS (o/w = 1/3) 


- 2 . 6 % 


. 2 % 






M HI t ( 


K ( % ERROR 
= +0.5 


FIGURE 14. ECCENTRIC CRACK (ISIDA'S PROBLEM) 
MODEL WITH TWO 8-NODE CRACKED 
ELEMENTS 



= 113 





FIGURE 15. 
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UNIFORM DISPLACEMENT, EQ. (81) 


YOUNG'S MODULUS 
E^ * 68.950 GPa 

E 2 = 30.683 GPa 

POISSON'S RATIO, 
v i *0.30 

X = 0 . 35 

STRESS INTENSITY FACTOR 
(REFERENCE 19) 

K 1 * 1.375S ] 

K 2 = 2.768S 2 

WHERE S IS REMOTE STRESS. 


CRACKED PLATE 





FIGURE 18. DISPLACEMENTS OF CRACKED FACES 
IN BI-MATERIAL PLATE 
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STRESS-INTENSITY FACTOR, K (MPa -/S) 



FINITE ELEMENT 
(NO RISER) 


68.95 MPa 


NO RISER 
LITERATURE 
(WILHELM) 
(REF 20) 


FINITE ELEMENT 
(WITH RISER) 


ASSUMING THAT FROM 
T * THIS POINT, THE RISER 
r AND PLATE CRACKS 

I RISER GROW AT THE SAME RATE 



CRACK LENGTH, a (mm) 


FIGURE 20. EFFECT OF RISER ON STRESS INTENSITY IN PLATE 




FIGURE 23. FINITE-ELEMENT REPRESENTATION OF THERMAL STRESS PROBLEM 


FIGURE A- 1. AXIAL ELEMENT 



FIGURE A-2. TRIANGULAR ELEMENT MEMBRANE 
DISPLACEMENTS 



FIGURE A-3. CONVENTION FOR DEGREES OF 
FREEDOM 





FIGURE A-5. DIRECTION VECTORS FOR FASTENER ELEMENT 
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FIGURE A -7. LOCAL COORDINATE SYSTEM FOR 
TEN-NODE CRACKED ELEMENT 


NASA- Langley, 1976 


CR-2698 



